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Abstract 

Creation of angular momentum in a relativistic electron-positron plasma is explored. It is shown 
that a chain of angular momentum carrying vortices is a robust asymptotic state sustained by the 
generalized nonlinear Schrodinger equation characteristic to the system. The results may suggest 
a possible electromagnetic origin of angular momenta when it is applied to the MeV epoch of the 
early Universe. 

PACS numbers: 52.60. +h, 52.40.Db 



1 



I. INTRODUCTION 



The problem of electromagnetic (EM) wave propagation and related phenomena in rela- 
tivistic plasmas has attracted considerable attention in the recent past. From the nonthermal 
emission of the high-energy radiation coming from a variety of compact astrophysical objects 
it has become possible to deduce the presence of a population of relativistic electrons in the 
plasma created in the dense radiation fields of those sources |JJ. The principal components 
of these plasmas could be either relativistic electrons and nonrelativistic ions (protons), or 
relativistic electron-positron (e-p) pairs. 

Relativistic e-p dominated plasmas may be created in a variety of astrophysical situations. 
The e-p plasmas are likely to be found in pulsar magnetospheres @ , in the bipolar outflows 
(jets) in Active Galactic Nuclei (AGN) ||, and at the center of our own Galaxy J|]. The 
presence of e-p plasma is also argued in the MeV epoch of the early Universe. In the standard 
cosmological model, temperatures in the MeV range (T ~ 10 10 K — lMeV) prevail up to times 
t = lsec after the Big Bang ||. In this epoch, the main constituent of the Universe is the 
relativistic e-p plasma in equilibrium with photons, neutrinos, and antineutrinos. 

Contemporary progress in the development of super strong laser pulses with intensities 
/ ~ 10 21_23 W/cm 2 has also made it possible to create relativistic plasmas in the laboratory 
by a host of experimental techniques |J. At the focus of an ultrastrong laser pulses, the 
electrons can acquire velocities close to the speed of light opening the possibility of simulating 
in the laboratory the conditions and phenomena that, generally, belong in the astrophysical 
realm [0. 

Elucidation of the electromagnetic wave dynamics in a relativistic plasmas will, perhaps, 
be an essential tool for understanding the radiation properties of astrophysical objects as 
well as of the media exposed to the field of superstrong laser radiation. Although the 
study of wave propagation in relativistic plasmas has been in vogue for some time, it is 
only in the recent years that the nonlinear dynamics of EM radiation in e-p dominated 
plasmas || has come into focus. The enhanced interest stems from two facts: 1) e-p plasmas 
seem to be essential constituents of the universe, and 2) under certain conditions, even an 
ultrarelativistic electron-proton plasma can behave akin to an e-p plasma ||. 

Recently, we firstly found dark soliton as well as vortex soliton solutions in e-p plasma 
|T0| . In Ref. [|ll|], it is also shown analytically that dark soliton is the natural nonlinear 
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coherent structure in unmagnetized cold e-p plasmas. However, it is conceivable that soliton 
solutions obtained in a one dimensional formulation will turn out to be unstable in higher 
dimensions, which may lead to the creation of vortex solitons. Since dark and vortex solitons 
are asymptotically nonvanishing, they have received much less attention than their localized 
cousins due to the generally accepted requirement that the fields be localized in physical 
system. However, in recent experiments studying laser field dynamics in different kinds of 
optical media, it was demonstrated that dark and vortex solitons can be readily created as 
superimpositions upon a localized field background Vortex soliton solutions are also 
found in imperfect Bose gas in the superfluids and they are extensively investigated and 



discussed |14). In this paper, we systematically investigate the instability of dark solitons 
and show that it can lead to the formation of vortex solitons. We speculate about interesting 
application of the vortex soliton in the early universe. 



In the recent paper ||10|| , we had developed an argument for the creation of domains of 
nonzero angular momentum in the MeV era of the early universe when it is supposed to be 
dominated by a plasma of e-p pairs. We first showed that in such a plasma the dynamics of a 
pulse of electromagnetic radiation, with a frequency much larger than the plasma frequency, 
is controlled by a generalized nonlinear Schrodinger equation (GNSE) with a defocusing 
nonlinearity. Then, borrowing a result of nonlinear optics [I2~| where the standard nonlinear 
Schrodinger equation (NSE) with a cubic nonlinearity has been investigated in great depth 
and detail, we conjectured that even the GNSE whose nonlinearity is similar in nature to 
that of NSE will allow vortex soliton solutions with an angular momentum that is conserved 
during propagation. The latter system also allows dark soliton solutions in one dimension 
which are known to be unstable to two-dimensional (2d) perturbations and eventually evolve 
to 2d vortex chains. 

In this paper we demonstrate that the arguments given in Ref. |10[ are actually borne 
out by the direct solutions of the derived GNSE. We begin by analytically showing that the 
Id dark solitons of the GNSE are, indeed, unstable to 2d perturbations. In the process, 
we derive the instability criterion which agrees with the numerical solution of the linearized 
system. 

The main part of the demonstration, however, comes from a numerical simulation of our 
GNSE. Starting from a broad variety of initial conditions, we find the eventual emergence 
of angular momentum carrying vortex solitons. That is, the stipulated solitonic structures 
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are readily accessible within the framework of this general equation with a nonlinearity 
more complicated than that of NSE. Thus we can state with much greater confidence that 
electromagnetism, operative in the MeV era, could easily be the primordial source of angular 
momentum associated with various structures of the observable universe. 

We again stress that the present results are quite general in the nonlinear dynamics in 
an electron positron plasma. 



II. SINGLE VORTEX SOLITON 



In the envelope approximation, a finite amplitude, circularly polarized electromagnetic 
(EM) pulse propagating in a relativistic e-p plasma obeys (for details, see [fT0[| ): 

2\ud T A L + c 2 VlA ± + - 2 dlA ± + ~)a ± = 0, (1) 

where A± is the slowly varying amplitude of the perpendicular (to the propagation direction) 
vector potential, V 2 = <9 2 + <9 2 , and £ = z — v g t is the "comoving" (with group velocity v g ) 
coordinate. Here to e = (47rno e e 2 /m e ) 1 / 2 is the electron Langmuir frequency, no = no e = uq p 
is the number density of the unperturbed background (subscript denotes the value at 
infinity), 7 is the Lorentz factor, and G = Ks(m e c 2 /T) / 'K2(m e c 2 /T) with K n , denoting the 
n-th order modified Bessel functions of the second kind, is the enthalpy density. 

In this approximation the continuity equation, determining the number density n, be- 
comes 

7 Gn [72Gg-G 2 -e 2 |AJ7(m e c 2 ) 2 ]V2' 1 ' 

and the system is closed by the adiabatic equation of state 

——^ — — exp(-m e c 2 G'/T) = const. (3) 

'yK 2 {m e c 2 /T) 

The low frequency motion of the plasma is driven by the ponderomotive pressure [ ~ 
(Pep) 2 ] °f tli e high-frequency EM field, and is independent of the sign of the particles' 
charge. It is perfectly natural to assume that the electron and the positron fluids have equal 
temperatures (T 0eiP = T ) in equilibrium so that their effective masses (G £tP = G) will also 
be equal. The radiation pressure will impart equal low-frequency momenta to both fluids 
allowing the possibility of overall density changes without producing charge separation. The 
charge neutrality conditions n e — n p — n, <fi — has been assumed by neglecting the small 



inequality of the charge due to baryon asymmetry. It is also evident that the symmetry 
between the two fluids keeps their temperatures always equal (T BiP = T) if they were equal 
initially In deriving (JI|)-(fJ), we have also assumed that the plasma is transparent (i.e. 
uj ^> uj e ), and that the longitudinal extent of the pulse is much shorter than its transverse 
dimensions (Ly <C L±). 

Defining the normalized variables 

^ uJ e el,. , iS 

t->t, -^=r^r, A ± ^A ± , (4) 



2ojG ' Cy/2G ' m e c 2 7 s G 

and noting that even with the assumption 3> Vl, the diffractive term can be the same 
order or even greater than the dispersive one for a highly transparent plasma (jg ^> 1), 
Eq. (0) converts to the following GNSE ( lg G > G): 

1 



id t A + -VlA + f(\A\ 2 )A = 0, (5) 



with the nonlinearity 



Equation (^|) admits a symmetric two-dimensional solitary wave solution. For stationary 
solitons, the ansatz 

A = A(r) exp(im9 - iXt), (7) 

reduces (||) to 



d 2 , irt/ IdA m 2 

ar z r ar r 

with the potential V(A) defined as 



> A + V 'W = --->- + -* A (8) 



V(A) = (A + 2) A 2 + 4Vl - A 2 - 4, (9) 

with the prime denoting the derivative with respect to its argument, and A representing the 
nonlinear frequency shift. 



A numerical analysis of the solution of this eigenvalue problem is given in Ref. ||10 
Because of the absence of ^-derivatives in @, any pulse-like localized structure is allowed in 
the propagation direction. In the direction perpendicular to the propagation, on the other 
hand, the amplitude was shown to approach a constant value. Since there does exist a small 



baryon asymmetry in the early Universe that acts on a longer characteristic length |15|, [16 
the extent of the constant amplitude region may be considered to be finite. 



Due to the single-valuedness of the vector potential, m must be an integer and A must 
vanish at the origin for non-zero values of m. The nonzero m solutions are particularly 
important because they carry the orbital angular momentum M 

(M) z = i Jdr x [r ± x (A*V ± A - c.c.)] z . (10) 

It is straightforward to show that the system (|5|) conserves angular momentum, and the 
expression (|l0|) is just the paraxial approximation for the orbital angular momentum, 
Me = J dr [r x (E x B)], of the EM field |fT7f| . The angular momentum carried by the 
vortex is M z = mN where the photon number N = J drj_ \A\ 2 is another conserved quan- 
tity; m is also known as the "topological charge." Strictly speaking, one must redefine the 



integrals of motion for non- vanishing boundary conditions |P-8| , but such a renormalization 
is not important here because of the fact that infinite-extent solution is just a formal ap- 
proximation. The presence of a small fraction of ions makes the physical solution decay at 



infinity M, 16 



III. DYNAMICS OF ANGULAR MOMENTA CREATION 



It is already suggested in Ref. [I(J that dark stripe solitons are unstable and break into 
vortex filaments. When the amplitude is small (\A\ <C 1), the nonlinearity reduces to the 
simple Kerr-type, and the vortex dynamics of NSE with such a nonlinearity has been studied 
in a variety of numerical calculations [|12|, [18], [19], |2(J. Our full nonlinearity is somewhat 



different and it has to be independently investigated. For this section, our goal is to show 
that the one- dimensional dark stripe solution is unstable against transverse perturbations 
in two dimensions. Such an instability causes the breakup of the stripe leading to a chain 
of vortex solitons with alternating polarity. 

A. Transverse instability 

The one-dimensional stationary solution for the GNSE (|5]) has the form g(x)e~ lXt (A > 0), 
where the real function g(x) satisfies 

*9 + ~§ + f(9 2 )9 = 0- (11) 



2dx 
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Since we are examining the linear stability of the solution with the least number of nodes, 
g is assumed to be an odd function with a single node at the origin. 
Let us perturb the one- dimensional solution by 

A = (g + u + w)e- iXt , (12) 

where the small perturbations u and v are real functions of x and y. 

Assuming sinusoidal behavior, u oc cos(ky — Qt) and v oc sin(ky — fit), the linearized 
GNSE could be reduced to the following eigenvalue problem: 

k 2 k 4 
Q 2 v = C 1 C v-—(C + C l )v + —v. (13) 

where 

1 8 2 



£i = X + \-^ + f(g 2 ) + 2g 2 f(g 2 ). (15) 



Remembering that Cog = and C\g x = 0, we can construct from (O) the relation 



14 i 2 7,2 

V 2 (g x | v) = -(g x | v) + (g x | dCtio) - -(g x \ 2C t v) - -(g x \ - 2g 2 f'v) (16) 

by multiplying g x on both sides and integrating with respect to x. Self-adjointness of Ci 
reduces ([TBI) to 

n2 _ k* k 2 (g x | - 2g 2 f'v) 

4 2 (g x \v) ' (17j 

provided (g x \ v) exists. A necessary condition for an exponential instability, then, is 

>0. (18) 



(gx 


- 2g 2 fv) 


{9x 


v) 



If it is true, then for sufficiently small \k\, we would have an instability, i.e., (0) is also a 
sufficient condition. Moreover, it is explicitly shown that uo — > o in the limit k — > unless 
(g x \v) diverges. 

We now make an approximate estimate of the ratio on the left hand side of (|T8|). Let us 
assume that the operators £ an d C-i are of order unity and Q 2 and k 2 are, in some sense, 
small. In this limit, we could estimate v by solving 

dC Q v = (19) 
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which has the solution v = g + q, where 



£o<? = 9x- 



(20) 



Notice that for v = g, the integrals occurring in QT7D are zero and this part will not contribute 
to the integral; only q will. From (|20|) , we can deduce the following: the first consequence 



(9x I g x ) = (9x | £o<?) = (9x | (£o - Ci)q) 

= (g x | - 2g 2 fq) = (g x | - 2g 2 fv) 



converts (O) into 



f? 2 



_ k 2 (g x | g x ) 
4 2 



(21) 



(22) 



To derive the second, we notice that 

(g x I v) = (g x | q) = (g x | C$ l g x ) (23) 

In principle, one can evaluate (p3|) in a rather straightforward way. But a very approximate 
estimate can be made by simply dropping d 2 in £ so that 

&c \ 



and 



a 1 = — 



k A k 2 



(gx I gx) 



2 (gx\g x /(X + f)) 



^ i 2 

ak . 

4 



In (25) we have overestimated the denominator. The eigenvalue is 



i7 = iak 



1 - 



4a 



1/2 



(24) 



(25) 



(26) 



implying a window < < 2a 1//2 in fc, where instability is possible . From the property 
that the one- dimensional solution is stable to one- dimensional perturbation, the growth rate 
Q approaches zero as the wave number k tends to zero. 

In order to accomplish the numerical analysis, we introduce a complex-valued function 
w = u + iv in (12). Then, the linearized equation looks 

1. 



id t w = Xw + ^Vlw + f(g 2 )w + 2g 2 f\g 2 ) Re(w), 

where the prime denotes derivative with respect to its argument. Putting dt 
d y = ik yields a linear eigenvalue equation for w: 

1 fd 2 w 



Qw = \w + 



2 I dx 2 



k 2 w +f(g 2 )w + 2g 2 f'(g 2 )Re(w). 



(27) 
-if} and 

(28) 



We have numerically solved (|28|) by the shooting method with boundary conditions w — > 
(x — > oo) and dw/dx — (x — 0). It is clear from fl28|) that the complex conjugate of the 
eigenfunction is also an eigenfunction. Namely, if w e is an eigenfunction of fl28|) corresponding 
to the eigenvalue i7 e , then w e is also an eigenfunction of ( |2"%| ) corresponding to the eigenvalue 
i? e , where the bar denoting the complex conjugate. 

For x — > +oo, g approaches a constant value go, which can be expressed as 



g ° = f-(—2)- (29) 
where the inhomogeneity due to the potential vanishes. In this region, x also becomes an 
ignorable direction, and we can assume w oc e~ KX for the point spectra (k > 0). Plugging it 
into our eigenvalue problem, we obtain the following "dispersion relation" which is applicable 
for large x, 



Of 



a + l^ 2 -k 2 ) + 2g 2 f'(g 2 ) 



(30) 



where the subscripts r and i denote the real and imaginary part, respectively, and we have 
used the relation f(gl) = —A. By solving it for k, we obtain 



K = yjk^-2Q t - 2g 2 f>(g 2 ) ± 2y/g$[f' (g 2 )] 2 - Q 2 , (31) 
where f'(g 2 ) is explicitly shown as 

/'(«) = -(1 - «)" 3/2 . (32) 
The asymptotic argument of w is also determined by 

With these informations, we have solved ( p8|) numerically. First, we fix an arbitrary value 
on the real part of w at x — 10. Given an expectant eigenvalue, the imaginary part of w 
is calculated from (|33"1) . Then, the derivative dw/dx is estimated by using (|3~lf) , and we can 
integrate (£8|) from x = 10tox = 0by the 4-th order Runge-Kutta formula. The data of g(x) 
is drawn from the numerical solution of (|TT]) . At x = 0, the boundary condition dw/dx = 
is checked. If this boundary condition is not satisfied, we guess the next eigenvalue by the 
Newton method and carry out the shooting again. 

The dispersion relation is shown in Fig. |I[ We can see a good qualitative agreement 
with the analytic evaluation with the fact that the growth rate begins from zero as k —>■ 0, 



experiences a maximum value with respect to k, and the mode is finally stabilized for a 
sufficiently large k. The eigenfunctions corresponding to k — 0.1, 0.5, and 0.7 for the 



the eigenfunction becomes wider and wider as k — > 0, which agrees qualitatively with the 
analytical estimate. Since we are performing the numerical shooting in the finite domain, it 
becomes difficult to estimate the correct eigenvalue in this regime. 

B. Nonlinear evolution 

We will now present the numerical simulation of the dynamics of angular momentum 
creation from one-dimensional dark stripe solitons by solving the full nonlinear system (|5])- 
(H). The calculation was carried out on a 200 x 200 spatial mesh placed in a calculation 
box of size L = 28. The boundary conditions d x A = (d y A = 0) are imposed at the edges 
x = ±L (y = ±L). Note that these boundary conditions numerically respect the constants 
of motion in the domain of integration. As an input we choose the one-dimensional x 
dependent stationary solution with A = 0.5. It is to be noted that the photon number and 
the Hamiltonian are conserved to order 10~ 5 during the calculation. Since the initial value 
of the integrated angular momentum is zero, it is found to always remain zero (within an 
error of order 10~ 15 ). 

Strictly speaking, the eigenfunction in the previous section (see Fig. |2|) contains cores 
of angular momenta, i.e. crossings of the real and imaginary zeros of the field due to the 
periodic form of the perturbation in the y direction. Thus, the growth of small amplitude 
perturbation itself does not exactly mean the 'creation' of angular momenta. In order to 
check the real creation of angular momenta from an exact zero everywhere, we have first 
carried out the calculation with an initial condition of the form 



with a = 0.9, (5 = 0.0875, and Ax = 0.1. This initial condition is so arranged that there is no 
crossing of the real and imaginary zeros. It is done by imposing a real- valued perturbation on 
a complex-valued stationary solution. The zero line of the imaginary part exactly coincides 
with the y-axis, while that of the real part deviates from the y-axis by A\. As the system 



parameter A = 0.5 are illustrated in Fig. 0. From Fig. ^ the mode with k = 0.5 gives 
the maximum growth rate Im(i?) ~ 0.137. As is seen from Fig. |2], the imaginary part of 




(34) 
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evolves, we observe the creation of crossings and the appearance of vortex solitons. As 
the initial conditions are changed (always starting form no crossings), the main qualitative 
result — the emergence of angular momentum carrying vortex soliton chains — is found to 
be fairly robust. 

Next we investigate the evolution of the system under transverse perturbations of the 
form 

A(x,y,0) = g(x)+A^coa(j^\, (35) 

with A = 0.1 and n — 6. The zero lines of the real and imaginary parts at t — 2 are shown 
in Fig. |3|. As we noted in Sec. O, the crossing points of two zero lines correspond to the 
vortex centers. We initially have twelve crossings suggesting twelve vortices. These twelve 
vortices, however, do not have well-formed solitonic structures since they are too close, and 
overlapping. 

As the system evolves, the vortices move and we see annihilations of pairs of opposite 
polarity. The first annihilation event was observed in the interval 28 ;$ t < 34; two of the 
vortices destroy one another near the center, and the other two seem to disappear near the 
edge; the annihilation of two pairs removes four from the original twelve at t — 2. 

The second annihilation event occurs during 40 < t < 48. As is depicted in Fig. |], the 
zero lines of the real and imaginary parts tend to separate around the center y ~ 0. Here 
two vortices approach the origin around 40 ;$ t ;$ 44, and then they are annihilated. After 
the annihilation, two inner vortices approach the origin. Finally the remaining six vortices 
tend to spread and align with equal inter- vortex spacing (see t = 72 in Fig. ^). 

The amplitude \A\ corresponding to this sequence is illustrated in Fig. |5]. At t — 48, 
the vortices do not quite look like vortex solitons, even though the number of vortices has 
been reduced to six. Since the distance between adjacent vortices is still small, they overlap 
and are not quite independent. The central hump around x ~ y ~ is the remnant of the 
second annihilation event. With time, the peaks tend to expand to the central region and 
are sufficiently apart to look and behave like vortex solitons. By t ~ 72, we observe the 
formation of six solitonic structures as predicted in Ref.|l0| (see t = 72 in Fig. ||). Notice 
that the vortices disappear only when annihilated by another with opposite polarity. 

After t •> 80, a propagating wave is clearly seen in the region y < — 10 and y > 10 where 
the field was originally flat (see Fig. |6]). This propagating field is the trace of the Cherenkov 



radiation which comes from the non-integrability of the system ||20|| . The radiation will 

11 



propagate away from vortices. In the finite calculation domain, the radiation will be reflected 
back at the boundary. 

When we extend the time evolution further, we observe the third annihilation event 
around t ~ 160. However, this annihilation may be an artifact of the finite size of the 
domain. When we carry out the simulation in a domain larger than (L = 28) but with the 
same initial conditions, the third annihilation event takes place at a later time while the times 
of the first and the second annihilation events remain unchanged. The third annihilation 
event is likely to be driven by the reflection of the Cherenkov radiation at the boundary. 
Equivalently, it may be concluded that six vortex solitons are stable in our domain and may 
stay forever. It is noted that the distance among vortices at the final stage approximately 
coincides with the inverse of the wave number with the maximum growth rate k m . 

IV. SUMMARY 

We have demonstrated the dynamics of angular momenta creation in a highly relativistic 
electron-positron plasmas subject to the passage of a strong pulse of electromagnetic fields. 
The system is governed by a generalized nonlinear Schrodinger equation with a defocusing 
inverse square root type nonlinearity. It turns out that the one-dimensional dark soliton 
stripe solutions of this equation, just like those of the standard nonlinear Schrodinger equa- 
tion, are unstable to transverse perturbations. By carrying out a linear analysis, we have 
found that there exists an instability window of transverse wave numbers for the system. 
By a numerical simulation of the fully nonlinear equation, we have shown that the trans- 
verse instability will yield, after a few annihilation events, a well-separated chain of vortex 
solitons with alternating, singly- charged polarity or topological charge (m = ±1). These 
singly-charged vortex solitons are topologically stable and do not disappear unless they col- 
lide with their compliments and annihilate. The number of the created vortex solitons seem 
to be determined by the inverse of the wave number k m with the maximum linear growth 
rate. For the box size L = 28, the six vortex soliton state is found to be robust. 

We have suggested a simple and plausible mechanism of angular momentum generation in 
the MeV epoch of the Universe. Electromagnetism, operating through the versatile substrate 
of the electron-positron plasma, seems to readily generate highly interesting, long-lived ob- 
jects which are capable of carrying large amounts of mass, energy and angular momentum. 
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Since an initial localization of mass, energy and angular momentum is precisely the seed 
that gravity needs for eventual structure-formation, electromagnetism may have provided a 
key element in the construction of the large-scale map of the observable Universe. 
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Figure Captions 

1. Dispersion relation for A = 0.2, 0.5, and I. 

2. Eigenfunctions for A = 0.5 corresponding to k = 0.1, 0.5, and 0.7. 

3. Zero lines of real and imaginary parts of the field at time t = 2. 

4. Crossing of zeros at time t = 40, 44, 48 and 72. 

5. Amplitude \A\ at time t = 48 and 72. 

6. Radiation appears after t > 80. 
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FIG. 1: Dispersion relation for A = 0.2, 0.5, and 1. 
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FIG. 3: Zero lines of real and imaginary parts of the field at time t = 2. 



t = 40 



t = 48 



t = 44 



t = 72 



FIG. 4: Crossing of zeros at time t = 40, 44, 48 and 72. 
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FIG. 5: Amplitude \A\ at time t = 48 and 72. 




FIG. 6: Radiation appears after t > 80. 
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